#Run stay-at-home order regressions using did package

rm(list=ls())

getwd()
setwd("C:/Users/g1stm01/Dropbox (Research)/stay_at_home_project")

#devtools::install_github("jonathandroth/staggered", upgrade="always")
#devtools::install_github("bcallaway11/did", upgrade="always")

library(tidyverse)
library(fixest)
library(broom)
library(staggered)
library(data.table)
library(did)
library(stringr)
library(haven)

set.seed(342021)

#Load data from Stata

reg_data <- read_dta("data/clean/statedata_foreventtime.dta")



##########Step 2: Use did packages to test out our results

reg_data_clean = reg_data %>%
  mutate(t=as.integer(date)) %>%
  mutate(daytreated=as.integer(effective_date)) %>%
  as.data.frame()

reg_data_clean = reg_data_clean %>% 
  mutate(t = replace_na(t, Inf)) %>% 
  arrange(state_fips, t)


#Loop over all outcomes

outcomes = c("sales_sector_percap_accommod", "sales_sector_percap_amazon", "sales_sector_percap_artsen",
             "sales_sector_percap_construc", "sales_sector_percap_educatio", "sales_sector_percap_finance",
             "sales_sector_percap_fooddel", "sales_sector_percap_general", "sales_sector_percap_grocery",
             "sales_sector_percap_healthc", "sales_sector_percap_hotelbo", "sales_sector_percap_informat",
             "sales_sector_percap_manufact", "sales_sector_percap_nonreta", "sales_sector_percap_otherse",
             "sales_sector_percap_overall", "sales_sector_percap_pharmaci", "sales_sector_percap_professi",
             "sales_sector_percap_realest", "sales_sector_percap_restaura", "sales_sector_percap_retail",
             "sales_sector_percap_retailt", "sales_sector_percap_transpor", "sales_sector_percap_utilitie",
             "sales_sector_percap_wholesal")
outcomes = c(outcomes,paste0("log",outcomes))


for (outcomevar in outcomes) {
  
  regdata_complete = reg_data_clean %>% select(state_fips,t,daytreated,state_fips,pop_estimate_2018,date,eval(outcomevar))
  regdata_complete = regdata_complete %>% filter(complete.cases(.))
  
  out = att_gt(yname = outcomevar,
               gname = "daytreated",
               idname = "state_fips",
               tname = "t",
               xformla = ~1,
               data = regdata_complete,
               weightsname = c("pop_estimate_2018"),
               allow_unbalanced_panel=TRUE,
               control_group=c("notyettreated"),
               est_method = "reg",
               clustervars = c("state_fips")
  )
  
  es <- aggte(out, type = "dynamic", min_e = -10, max_e = 10)
  summary(es)
  
  
  estimates <- es[["att.egt"]]
  se <- es[["se.egt"]]
  relative_time <- es[["egt"]]
  
  graph_data <- data.frame(relative_time, estimates, se)
  
  overall_att <- round(es[["overall.att"]], digits = 3)
  overall_se <- round(es[["overall.se"]], digits = 3)
  
  
  base_graph <- graph_data %>% 
    mutate(ymin_ptwise = estimates + 1.96*se,
           ymax_ptwise = estimates - 1.96*se) %>% 
    ggplot(aes(x=relative_time, y =estimates)) +
    geom_point() + 
    geom_ribbon(aes(ymin=ymin_ptwise, ymax=ymax_ptwise), alpha=0.2) +
    geom_hline(yintercept =0, linetype = "longdash", alpha = 0.15) +
    geom_vline(xintercept = 0, linetype = "longdash", alpha = 0.15) +
    scale_x_continuous(breaks = seq(-10, 10, by = 1), minor_breaks = seq(-10, 10, by = 1)) +
    xlab("Days relative to order") +
    annotate(geom = 'text', label = paste0("Overall ATT = ", overall_att), x = -Inf, y = Inf, hjust = 0, vjust = 1) +
    annotate(geom = 'text', label = paste0("Overall SE =  ", overall_se), x = -Inf, y = Inf+1, hjust = 0, vjust = 2.5) +
    theme_minimal() +
    theme(axis.text.x=element_text(size=rel(1.2)),
          axis.text.y=element_text(size=rel(1.2)))
  
  plot(base_graph)
  
  setwd("C:/Users/g1stm01/Dropbox (Research)/stay_at_home_project/figures/did")
  ggsave(filename = paste0("sm_",outcomevar,".png"), dpi = 600)
  
}



